function [] = ROMS_tidegauge_extractor(indirlist,mlocations,outpath)
% ROMS_tidegauge_extractor.m
%
% Essentially a "driver" script for Z_get_tidegauge.m to extract limited
% information from ROMS history files near specified tide gauge locations.
% Very similar to the "mooring extractors" except more limited information
% is extracted since actual tide gauges only record limited information.
% Based on previous code by S. Giddings, but updated to avoid looping
% through the locations one-by-one by making use of the vector lat, lon
% ability of Z_get_tidegauge.m.
%
% NOTE:  The codes extracts output at the nearest wet rho-grid point, not
% the exact tide gauge position, because the actual tide gauge position may
% be covered by the model land mask. However, Z_get_tidegauge.m is still
% calculating weights for interpolation even though we don't need to
% interpolate...
%
% REQUIRES:  tidegaugeLocations.mat - which is a .mat file with lon, lat
% positions for the named tide gauge locations.  See Z_tidegaugeLocations.m
% and make_tidegaugeLocations.m to update the .mat file.
%
%
% TO RUN ON SKUA OR FJORD IN THE BACKGROUND:  in a terminal window go to
% the directory where the function lives and type
%
%     matlab < ROMS_tidegauge_extractor.m > & TG_extract_log &
%
% the typical extraction for one year's worth of hourly history files takes
% roughly 45 minutes on 'fjord' for all tide gauge locations.
%
%
% R. McCabe 04/2014


%% User input to choose the run to work on
disp(' ');
disp('********** Tide Gauge Extractor **********');
disp(' ');

% Set list of input directories to loop through
if ~exist('indirlist')
    % indirlist = {'put your path here, or multiple paths here!'};
    % indirlist = {'/pmraid4/sarahgid/runs/ptx_highT_2_2005/OUT','/pmraid4/sarahgid/runs/ptx_highT_2_2006/OUT','/pmraid4/sarahgid/runs/ptx_highT_2_2007/OUT'};
    % NOTE: leave off the '/' after 'OUT' in indirlist to get the correct basename below
    indirlist = {'/pmraid1/siedlesa/T40final_base_2005_noramp/OUT'};
end

% Set list of tide gague locations to extract
if ~exist('mlocations')
    
    mlocations = {'Charleston'
         'South Beach'
         'Garibaldi'
         'Astoria'
         'Toke Point'
         'Westport'
         'La Push'
         'Bamfield'
         'Tofino'
         'Neah Bay'
         'Port Angeles'
         'Victoria'
         'Port Townsend'
         'Friday Harbor'
         'Patricia Bay'
         'Cherry Point'
         'Vancouver'
         'Point Atkinson'
         'Seattle'
         'Tacoma'};
end

% Set the output directory for the extracted mooring records
if ~exist('outpath','var')
    % outpath = 'put your desired output path here!';
    outpath = '/data1/rmccabe/runs/PostProcess/ROMS_tides/';
end
if ~exist(outpath,'dir')
    mkdir(outpath)
end


%% Set global paths and directories
addpath('../../alpha/');
Tdir = toolstart;


%% Load tide gauge names and locations
load([Tdir.tools 'validation_tools/tide_gauge/tidegaugeLocations.mat']); % tgauge
% Alternatively see Z_tidegaugeLocations.m and make_tidegaugeLocations.m to
% change and resave the file tidegaugeLocations.mat


%% Loop through input directories and make the tide gauge extractions
tic
for inum = 1:length(indirlist)
    indir = indirlist{inum};
    
    % choose which run to use, and set the basename
    % assumes that "indir" is something like:
    % /Users/PM/Documents/Salish/runs/ptx_med_2005_1/OUT
    % then the lines below find the string right before "OUT", which is the
    % basename, the main identifier of the run
    % NOTE: leave off the '/' after 'OUT' in indir to get the correct basename
    ind = strfind(indir,'/');
    basename = indir(ind(end-1)+1:ind(end)-1);
    disp(' '); disp(['basename = ',basename]);
    
    % Get the run definition
    runi = roms_createRunDef('my_run',indir);
    % Get the year of the run
    yeari = datestr(runi.his.nctime(1),'yyyy');
    yeari = str2double(yeari);
    
    % Initial message
    disp(' '); disp(['Extracting tide gauges:']); disp(' ');
    
    % Preallocate NaN vectors for mooring lons,lats
    lon_moor = NaN.*ones(length(mlocations),1);
    lat_moor = NaN.*ones(length(mlocations),1);
    
    % Loop through locations to extract to build lon,lat vectors
    for mnum = 1:length(mlocations)
        mloc = mlocations{mnum};
        
        % Identify the current tide gauge location
        gi = find(strcmp(cat(2,{tgauge.station_name}),mloc));
        lon_moor(mnum) = tgauge(gi).lon;
        lat_moor(mnum) = tgauge(gi).lat;
        clear gi
        
        % List the tide gauges and locations for confirmation (not actually
        % extracting yet...)
        disp(['  Tide gauge ',mloc]);
        disp(['      Lon = ',num2str(lon_moor(mnum)),...
            '      Lat = ',num2str(lat_moor(mnum))]);
        disp(['    Locating nearest wet grid point...']);
        
        % Find the nearest point on the rho-grid, that is not on the mask!
        % NOTE: u,v-grid items are retained for completeness, but commented
        % out below. See roms_loadGrid.m for exact variable definitions.
        llon = runi.grid.lon; llat = runi.grid.lat;
        %llonu = runi.grid.lonu; llatu = runi.grid.latu;
        %llonv = runi.grid.lonv; llatv = runi.grid.latv;
        Llon = llon(:); Llat = llat(:);
        %Llonu = llonu(:); Llatu = llatu(:);
        %Llonv = llonv(:); Llatv = llatv(:);
        mask = runi.grid.mask;
        %masku = runi.grid.masku;
        %maskv = runi.grid.maskv;
        %
        % Just look for points where the mask == 1
        Llon(mask(:)==0) = []; Llat(mask(:)==0) = [];
        %Llonu(masku(:)==0) = []; Llatu(masku(:)==0) = [];
        %Llonv(maskv(:)==0) = []; Llatv(maskv(:)==0) = [];
        %
        % Search for the nearest lon/lat on the mask
        k = dsearchn([Llon,Llat],[lon_moor(mnum),lat_moor(mnum)]);
        %ku = dsearchn([Llonu,Llatu],[lon_moor(mnum),lat_moor(mnum)]);
        %kv = dsearchn([Llonv,Llatv],[lon_moor(mnum),lat_moor(mnum)]);
        %
        % Then find that same lon/lat in the full lon/lat matrix
        ixr = dsearchn(llon(1,:)',Llon(k));
        jyr = dsearchn(llat(:,1),Llat(k));            
        %ixu = dsearchn(llonu(1,:)',Llonu(ku));
        %jyu = dsearchn(llatu(:,1),Llatu(ku));
        %ixv = dsearchn(llonv(1,:)',Llonv(kv));
        %jyv = dsearchn(llatv(:,1),Llatv(kv));
        clear llon llat Llon Llat mask k
        %clear llonu llatu llonv llatv Llonu Llatu Llonv Llatv
        %clear ku kv masku maskv ixu jyu ixv jyv
        
        % Calculate approximate distance to nearest wet rho-grid point
        [mdist,~] = sw_dist([runi.grid.lat(jyr,ixr) lat_moor(mnum)],...
            [runi.grid.lon(jyr,ixr) lon_moor(mnum)],'km');
        
        % Update lon_moor lat_moor to nearest wet rho-grid point
        lon_moor(mnum) = runi.grid.lon(jyr,ixr);
        lat_moor(mnum) = runi.grid.lat(jyr,ixr);
        clear jyr ixr
        
        % Tell everyone what we're really doing!
        disp(['    Will extract output at nearby rho-grid point ',...
            num2str(round(1000*mdist)/1000),' km away:']);
        disp(['      Lon = ',num2str(lon_moor(mnum)),...
              '      Lat = ',num2str(lat_moor(mnum))]); disp(' ');
        clear mdist
    end
    clear mnum
    
    % Extract the tide gauge records...
    % For testing, you may only want info from say XX output files...
    % see line 120 in Z_get_tidegauge.m to change NT = XX
    [TG] = Z_get_tidegauge(runi,lon_moor,lat_moor);

    
    % Pack the extractions into individual mooring .mat files and save...
    for mnum = 1:size(TG,2) % should be same as 1:length(mlocations)
        mloc = mlocations{mnum};
        
        % Pack the results in a structure
        TG(mnum).run = runi;
        TG(mnum).basename = basename;
        % Rename for backward consistency
        ssh_moor = TG(mnum);
        ssh_moor.name = mloc;
        % Remove spaces for file name
        mloc_nospace = mloc;
        mloc_nospace(mloc_nospace == ' ') = [];
        % Save the output tide gauge file
        save([outpath,basename,'_',mloc_nospace,'.mat'],'ssh_moor');
        
    end
    clear mnum
    
end %end of indir loop
elapsed_time = toc;
disp(' **FINISHED**');
